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ABSTRACT 

We present a ‘two-fluid’ implementation of dust in smoothed particle hydrodynamics (SPH) 
in the test particle limit. The scheme is able to handle both short and long stopping times and 
reproduces the short friction time limit, which is not properly handled in other implementa¬ 
tions. We apply novel tests to verify its accuracy and limitations, including multi-dimensional 
tests that have not been previously applied to the drag-coupled dust problem and which are 
particularly relevant to self-gravitating protoplanetary discs. Our tests demonstrate several key 
requirements for accurate simulations of gas-dust mixtures. Firstly, in standard SPH particle 
jitter can degrade the dust solution, even when the gas density is well reproduced. The use 
of integral gradients, a Wendland kernel and a large number of neighbours can control this, 
albeit at a greater computational cost. Secondly, when it is necessary to limit the artiflcial vis¬ 
cosity we recommend using the |Cullen & Dehn^ ( |2010| ) switch, since the alternative, using 
a ~ 0.1, can generate a large velocity noise up to crj < O.Scg in the dust particles. Thirdly, we 
And that an accurate dust density estimate requires > 400 neighbours, since, unlike the gas, 
the dust particles do not feel regularization forces. This density noise applies to all particle- 
based two-fluid implementations of dust, irrespective of the hydro solver and could lead to 
numerically induced fragmentation. Although our tests show accurate dusty gas simulations 
are possible, care must be taken to minimize the contribution from numerical noise. 


Key words: hydrodynamics - methods: numerical - protoplanetary discs - planets and satel¬ 
lites: formation - dust, extinction 


1 INTRODUCTION 

The growth of planets out of km-sized planetesimals remains an 
attractive scenario for both terrestrial planet formation and, in the 
core accretion model, also as a necessary step in the creation of gas 
giant planets. Nevertheless - regardless of whether these planetesi¬ 
mals form by direct gravitational collapse or agglomerative growth 
- there are a number of problems associated with the formation and 
retention of planetesimals in protoplanetary discs (for a recent re¬ 
view of the subject see |Johansen et al.| ( |2014| )). 

For example, although there is observational evidence at mm 
and cm wavelengths for rather rapid grain growth ( [Wilner et al.| 
|2005[ |Rodmann et al.|[T006[ |Ricci et aL]|2010| ), agglomerative 
growth beyond ~ cm-sizes is inefficient, since collisions tend to 
result in fragmentation or bouncing rather than sticking and grain 
growth. Particles in this size range are said to be ‘critically cou¬ 
pled’ to the disc gas, having values of the Stokes number (St: the 
ratio of the timescale for gas drag to the dynamical time) of around 
unity. Drag causes such particles to lose angular momentum to the 
gas (which orbits at slightly sub-Keplerian velocity on account of 

* E-mail: rab200@ast.cam.ac.uk 


being partially supported by outwardly directed pressure gradients) 
and to drift towards the star on time-scales of 100 to 1000 orbits 
( |Weidenschilling|1977| l. Once the drifting particles reach the snow¬ 
line they will sublimate, enhancing the local density in vapour form 
( |Cuzzi & Zahnle|2004|l, removing most of the solids from the outer 
disc within 1 Myr ( [Brauer, Henning & Dullemond|2008] ). 

The loss of solids by radial drift can potentially be overcome 
if pressure does not increase monotonically with decreasing radius 
in the disc and if instead some process can create and sustain lo¬ 
cal pressure maxima. In this case solid material can be trapped at 
the maximum (since regions of outwardly increasing pressure drive 
outward migration of solids by the same argument). For example, 
the pressure maxima at the outer edge of gaps formed by massive 
planets have been proposed as an effective dust trap ( |Lyra et ^ 
|2009| ), and this has been invoked to explain the asymmetric dust 
around Oph IRS 48 ( [van der Marel et al.|20I3| l. During the early 
phase of protoplanetary disc evolution, while the disc is still mas¬ 
sive enough that gas self-gravity is important, spiral structures may 
likewise slow or halt the radial drift of dust grains, concentrating 
dust in the pressure maxima associated with spiral arms ( |Rice et af] 
|2004 1 [2006) [Gibbons, Rice & Mamatsashvili 201^ . If this focus¬ 
ing of dust into the spiral arms can raise its density by two orders 
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of magnitude then it becomes the dominant mass component lo¬ 
cally and direct gravitational collapse in the dust layer - creating 
km scale planetesimals - may ensue. 

If planetesimals do form in self-gravitating discs, then the spi¬ 
ral structure will drive eccentricities to e > 0.1. The high velocity 
dispersion (and hence suppression of gravitational focussing) ren¬ 
ders collisions negligible and hence planetesimals should survive 
the self-gravitating phase { |Walmswell, Clarke & Cossins||2013| l. 
Planetesimal formation is most likely to occur at radii greater than 
10 AU , where the relatively short cooling time is linked to a larger 
amplitude of spiral structure ( [Clarke & Lodato||2009| l. The direct 
conversion of grains that are critically coupled (with size 1 - 10 cm 
in this region of a self-gravitating disc) into planetesimals offers 
the attractive prospect of by-passing the problems of agglomera- 
tive growth over intermediate size scales. 

Since the mechanism requires a significant fraction of the dust 
to be in large grains that have a Stokes number, St 1, the success 
of this model depends on their survival during the self-gravitating 
phase. Their survival is sensitive to relative velocity of collisions 
between particles, since collisions with velocities gre ater than 1 


10 ms ^ result in fragmentation ( [Guttler et al.|2010 i. Since 


Rice 


jet al.j ( |2004| ) found velocity dispersions of order the sound speed, 
which is approximately 500 ms“^ at 30 AU, collisions may de¬ 
stroy the grains before a sufficient density has been built up for 
collapse to occur. 

The results of competing grain growth and fragmentation 
processes require detailed knowledge of the velocity distributions 
of the dust grains during the self-gravitating phase, and this re¬ 
mains uncertain. In particular, early simulations of self-gravitating 
discs were studied using Smoothed Particle Hydrodynamics (SPH), 
which had problems related to the need to reduce artificial viscosity 
in order to avoid a situation where the heating by the gravitational 
instability is exceeded by the action of artificial viscosity on Ke- 
plerian shear ( jLodato & Rice||2004[ [Rice et aL]|2004| [2006 j ). The 
amount of artificial viscosity in SPH is typically controlled by two 
parameters, a and (3, which govern artificial viscosity terms that 
are linear and quadratic in the relative velocity, respectivel 3 Q To 
reduce the viscosity often a was reduced by an order of magnitude 
from the typical values of a — 1 and sometimes /3 was also re¬ 
duced ( jLodato & Rice|200l^ [Rice et al.|2004||2006[[Meru & Bat^ 
[201 Ij |2012| ). Subsequently it has been shown that such a choice 
fails to generate enough entropy in shocks and results in consid¬ 
erable noise in the gas velocities ( Meru & Bate|201l] 2012 


jet al.||20141 see also [Lattanzio et al.||1986[|Lornbardi et al. 


Rice 


\m9 


and Thacker et al. H2000) who suggest a > 0.7, even for weak 

shocks). The velocity noise is likely to be transmitted to the dust 
component via drag forces and this means it is not entirely clear 
that the velocity dispersion measured in a simulation is physical. 
Recent shearing-box simulations go some way towards addressing 
the issue ( [Gibbons, Rice & Mamatsashvili|2012|[Gibbons, Mamat-j 
[sashvili & Rice[[2014[ t and also find a velocity dispersion of dust 
particles over the entire box of u ~ c^. However, the collision 
velocities are sensitive to the correlation structure in the velocity 
distribution, a quantity that has not been derived from these simu¬ 
lations. 

With recent developments in computational techniques, accu¬ 
rate simulations of the dynamics of gas-dust mixtures are within 


^ For a detailed description of artificial viscosity used in SPH, see [Price [ 
( [20T2) ; although note that the code used in this work (gadget-2) uses a 
slightly different form jSpringel[2005| . 


reach ( Johansen et al. [[20(371 [Laibe & Price[[201^ [2014[ [Loren- [ 

[Aguilar & Bate|20T4] ). However, while there have been significant 

efforts to verify codes involving the simulations of turbulence and 
the streaming instability (e.g. jBai & Stone|2010[ ), there have been 
few ways to verify that codes capture the velocity distribution of 
solids in the context of self-gravitating discs. 

In this work we conduct a series of tests that can be used to 
assess code performance in this respect. While some of these (e.g. 
the DUSTYBOX and DUSTYSHOCK tests described in lLaibe & Pricel 
( |2011[ )) have previously been used in this way, we also introduce 
two tests that have not previously been examined with respect to 
the modelling of coupled dust. These two tests involve shocks in 
2D: the implosion test of jHui, Li & Li[ ( [T9^ and the modelling 
of a rigidly rotating spiral mode in an isothermal gas disc ( [Roberts [ 
[1969 1 [Shu, Milione & Roberts|1973| l. The latter in particular was 
chosen because it bears some resemblance to the problem of inter¬ 
est (since it involves centrifugally supported gas interacting with a 
pattern of spiral shocks) but - unlike the ^^//-gravitating disc prob¬ 
lem - is amenable to an analytic solution for the gas that is steady in 
the co-rotating frame. For each test we quantify the velocity struc¬ 
ture in the dust and how this is affected by a range of numerical 
choices. Finally, we summarise our conclusions in terms of a set 
of recommendations and requirements for the simulation of dust in 
self-gravitating protoplanetary discs. 


2 GAS - DUST DYNAMICS 

The dynamics of solids in protoplanetary discs is governed by the 
interplay of drag forces with Keplerian motion. The force on a solid 
particle of mass rrid is given by 


( 1 ) 


( 2 ) 


Fd = ^ -Ksi^d - Vg) + mdad,ext, 

and the corresponding equation for the gas is 
dv^ _ VP pd j. . . 

j , —-1- Ks^Vd — Vg,) + ag, ext, 

dt pg rridPg 

where and Vd are the local gas dust and gas velocities, pg, and 
Pd are their respective densities, P is the gas pressure and a^,ext 
and ad,ext are additional accelerations. These have been made dis¬ 
tinct for the gas and dust to allow for forces that affect both phases, 
such as gravity and those that only affect a single phase, like vis¬ 
cosity. In these equations ^ refers to the full Lagrangian derivative 
in the frame moving with the gas or dust. 

The drag coefficient. Kg, depends on the Reynolds number, 
gas density and temperature, which can vary significantly over the 
different gas and grain properties found in astrophysical problems. 
The stopping time, tg, can be defined as the time over which the 
velocity difference Av = ^d — decays, = — Av/ts, 

giving 


tg = 


rridPg 


( 3 ) 


Kg{pg + pd) 

For negligible dust densities, pdjPg ^ 1, this is equal to the stop¬ 
ping time for a single dust particle, rud/Kg. In protoplanetary discs 
an approximation for the drag coefficient by [ WhippT^ ( [ 1 972[ ) is gen¬ 
erally used. For small grains, where the grain size, a, is less than 
the mean free path. A, the drag force can be calculated using the 
Epstein approximation 


_ 47r 2 

Kg — ^ PgOj Vg, 


(4) 


where Vs = ^ySksT/nn, T is the gas temperature, fcs is the 
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Gas and dust in SPH 3 


Boltzmann constant and /x is the mean-molecular weight. For a > 
9 A/4, the drag force can be calculated in the Stokes regime, 


Ks = -CD7Ta^pg\vd - v^|. 


The coefficient Cd can be approximated by 


( 24/Re, 

Re < 1, 

Cd = 1 24/Re°®, 

1 < Re < 800, 

[ 0.44, 

Re > 800, 


(5) 


( 6 ) 


We take the approach of interpolating the gas properties to the loca¬ 
tion of the dust particles, since it is less sensitive to the underlying 
gas particle distribution. 

The gas properties are interpolated to the location of each dust 
particle using a kernel sum. For a property A, where each gas par¬ 
ticle has the value Ai, its value at the location of a dust particle Ad 
can be evaluated via 

= V^-^VF(|rj-rD|,ftD), (7) 

, Pj 


where. Re = 2apg\\' — | /zx is the Reynold’s number of the fiuid 

and u is the molecular viscosity of the gas. Typically in astrophys- 
ical applications pdjPg ^ 1. In these circumstances the effect of 
drag forces on the gas can be neglected and the dust particles can be 
considered as test particles, with pd — 0. While there are interest¬ 
ing astrophysical problems in which non-negligible dust density is 
required, such as the streaming instability ( [Johansen et al.|2007] l, 
the test-particle limit has the advantage of allowing comparison 
between the simulations and analytical solutions. Furthermore, it 
allows a more direct comparison with work on self-gravitating pro¬ 
toplanetary discs, in which the test particle limit was used (c.f. jRicej 
jet al.|2004||2006[ [Gibbons, Rice & Mamatsashvili|2012| l. 


2.1 Implementation 


A number of different approaches to including the mutual dynamics 
of gas and solids in SPH have been proposed, which can be broadly 
separated into single-fiuid or two-fluid approaches (Rice et si] 


[2004[[Laibe & Price|2012[ [20141 [Loren- Aguilar & Bate|2014[ ). In 

two-fluid approaches the gas and solids are integrated by including 
separate particle types for the gas and dust grains, while the single¬ 
fluid approach considers the dust fraction e = pdl{pg + Pd) and 
relative velocity Av = Vd — Vg to be properties of the SPH par¬ 
ticles and solves equations for these quantities on an unstructured 
mesh made up of the SPH particles. 

In the case of high dust to gas ratios, where momentum feed¬ 
back on the gas from the dust is important [Laibe & Price[ p014[ ) 
have shown that it is essential that all length scales in the problem 
are resolved for both the gas and dust treatments, which places a 
requirement on traditional two-fluid approaches that the smooth¬ 
ing length, h < tgCg, where Cg is the sound-speed. Single fluid 
approaches appear to avoid this difficulty, but they do so at the 
cost of the full phase-space information, since only one velocity 
each of the gas and dust components can be known at each point 
in space. Conversely, the two-fluid approach is able to capture full 
phase space information. In the test particle limit the motions of the 
dust particles are independent of each other and thus their motions 
can be calculated accurately as long as the relevant scales in the 
gas are resolved. Thus the criterion h <tsCs does not apply unless 
the gas varies on these scales. [Loren- Aguilar & Bat^ ( [2014[ ) have 
shown this idea holds in general for low dust-to-gas ratios. 

Since the full phase-space is essential for evaluating the col¬ 
lisions between grains and the growth of planetesimals, we take 
a two-fluid approach, which we implement into the SPH code 
GAD GET- 2 ( [Springel[|2005[ l. The two-fluid method has been im¬ 
plemented either by considering pair-wise forces between gas and 


dust particles (Monaghan & Kocharyan[ 1995 [[Laibe & Price[2012| 

[Loren- Aguilar & Bate 

2014J, or by interpolating the gas properties 

to location of the dust ( 

Rice et al.[2004[[2006|). When back-reaction 

is included pair-wise 

forces are preferable, since they explicitly 


conserve angular momentum. However, for the test-particle limit 
we find a similar performance for both methods (see Appendix [A|. 


where r^, rrij, pj and Aj are the position, mass and density and 
quantity to be interpolated evaluated at the location of the jth gas 
particle, and W is the smoothing kernel. The sum is over the gas 
particles neighbouring the dust particles and smoothing length, ho, 
is set in the same way as for the gas particles, by solving 

iVNGB = -ro|,/iD), (8) 

3 

for Hd Sit the location of the dust particle td - The factor ah^ is the 
volume of the kernel, in three dimensions zx = 3 and a = 47r/3; 
and Angb is the desired number of neighbours. The smoothing 
length and interpolated gas properties are evaluated for all active 
dust particles after the gas density has been calculated. As a default 
choice for the kernel, W{r, /i), we use the standard cubic-spline 
kernel and choose Angb = 50. In some tests we also consider 
Wendland kernels. We note this is contrary to recommendation of 
[Laibe & Price [ ( [2012[ ), but as shown by the tests in the following 
section, this choice generates the expected forces. 

2.2 Time Integration 

Time integration for gas particles in GADGET-2 uses a Kick-Drift- 
Kick leapfrog algorithm, in which the gas velocities and positions 
are updated alternatively. The positions, and velocities, at 
time t are advanced to time t St vm 


/n n , 

V = V + -a , 

(9) 

xu+i = x” + 

(10) 

^n+l ^ 

(11) 


where the total accelerations are evaluated at positions 
The leap-frog scheme is a symplectic, time-reversible integrator for 
Hamiltonian mechanics, which has the advantage that for closed 
orbits the energy error can neither grow or decrease, resulting in 
excellent long term stability even when the truncation error of a 
single step is relatively high. 

For dust particles with St ^ 1, which follow weakly per¬ 
turbed Keplerian orbits, the leap-frog integrator is a natural choice. 
However, for St 1 explicit time integration schemes require very 
small time steps. Since in these cases tg is much shorter than any 
other time-scale in the problem, this suggests the use of implicit, 
or semi-implicit approaches that allow At ^ tg. [Loren- Aguilar &| 
[Bate[ ( |20T^ extend this idea by using the analytical solution in an 
operator-split approach. They suggest initially updating the veloc¬ 
ities according to non-drag forces and subsequently applying the 
drag. The first step is 

v2 = Vcz + ad At, 

VP 

V* = + dig At - 

where from here we drop the explicit ext on ad,ext and a^,ext- The 
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drag forces are applied making use of the analytical solution for rel¬ 
ative velocity, Av, which can be derived in the absence of external 
forces. For linear drag laws this gives 

Av(t + At) = Av(t) exp(—At/ts). (12) 

Writing ^(t) = 1 — exp(—t/ts), this generates the final dust ve¬ 
locity via 

Vd(t + At) = v3-^(At)(v3 -V*), (13) 

which is valid in the test particle limit. At ^ ts leads to 
V(i(t + At) ^ V*. This property means that although the scheme 
is stable for all time-steps, it does not produce the correct termi¬ 
nal velocity unless At <C ts. To see this consider a static atmo¬ 
sphere in which the external forces include only constant gravity, 
ag = ad = g. For a static atmosphere both gas velocity and accel¬ 
eration must be zero, with pressure balancing gravity, VP/p = g. 
For large time-steps ^ ^ 1 and Vd(t + At) ^ 0, which is not the 
analytical solution, ^ t^VP/p. 

Instead, we approximate Wg as constant in both time and space 
during the time-step, giving rise to the dust kick 

Vd(t + At) = Vd exp(-At/ts) + (adts + v^)^(At), (14) 


which is exact for static atmospheres and uniform flow for all 
At /ts and reproduces the short friction time limit. This approxi¬ 
mation is as good as the approximation that and au are con¬ 
stant. Only when both the pressure gradients and external forces are 
zero do the analytical solution and operator split time-steps produce 
equivalent kicks. 

In some cases it may be beneficial to take into account the 
effect of motion of the gas on the update for Av (since uniform 
gas motion is already taken into account in equation [^. The full 
equation of motion for Av in a frame moving with the dust is 


dAv 

dt 


Av 

ts 


/ VP, 

+ (ad - + ——) 

Pg 


+ Av Vv^, 


(15) 


where the last term takes into account the relative advection be¬ 
tween the dust and gas (see equation 27 of |Laibe & Price] p012| )). 
While the advection could be calculatec0 neglecting it and approx¬ 
imating the forces as constant during a time-step gives 


VP 

Vd(t + At) = Vd exp(-At/ts) + Vg^{At) + (a^- )At 

P 


+ ad — ao + 


VP \ 
pg ) 


p^(At). 


(16) 


The utility of this approximation can be seen by looking at the lim¬ 
its ts ^ 0 and ts oo. In the first limit, ^ ^ 1 and the dust ve¬ 
locity reduces to the short friction time limit, but using the updated 
gas velocity. For the ts — 0, the dust velocity becomes exactly the 
updated gas velocity. In the second limit, ^ At/ts and the dust 
kick reduces to the explicit update for the dust. Since the approx¬ 
imation of constant ad, a^ and VP/p better refiects assumptions 
made for the integration of the equations of motion for the gas, we 
recommend the use of equation However, in practice we have 
used equation[^in the results presented. Our tests have shown that 


^ While there may be situations in which the relative advection term is 
the dominant source of error in equation there are many situations in 
which equation [Tb] represents an improvement. In the cases where the rel¬ 
ative advection term is important, it will be adequately controlled by the 
Courant-like condition because the velocity should not vary significantly 
on scales of order of the smoothing length in resolved flow. The exception 
to this is at shocks, where all schemes are limited to first order anyway. 


the difference in the results presented is small, because in general 
the SPH forces and interpolation dominate the error. 

To demonstrate the accuracy of the time-stepping schemes, we 
show the L 2 error norm for the velocity in Fig. for two tests in 
which we have prescribed a total gas acceleration a = VP/p + a^. 
In the first we used a ID forced oscillation, a = Vu cos{ijjt) and 
we integrated both the gas and dust using the leap-frog schemes. 
For the second test we used a travelling wave a = Vusinjut — 
kx), but integrated only the dust using the leap-frog schem^ The 
two tests are chosen to highlight the effect of neglecting the relative 
advection term, which is identically zero in the first example. For 
the parameters we have used V = 10~^, k = w = 1, and a 
stopping time of ts = 0.05. For the L 2 error norm, we use 


L 2 = 


J2iiri - vo{xi,t)) 


2 11/2 


Y,vo{xi,t) 


(17) 


where, for the forced oscillation we used the analytical solution 
for uo, and a high accuracy numerical solution for the travelling 
wave. For the range of number steps shown, both the regimes At <C 
ts and At ^ ts are tested, which corresponds to A ^20 and 
A 20 for these parameters. We see that in both cases equation 
[^produces a first order scheme, while equationis a second 
order scheme. 

If an explicit update for the dust velocity had been used, then 
for A < 50 the integration would have lead to significant error. 
However, as the stopping time decreases the power of our scheme 
becomes apparent, since we find essentially the same error regard¬ 
less of stopping time with our scheme, and an accuracy of 0.1 per 
cent can be achieved with 100 steps. Many more steps would be 
needed for tightly coupled particles (ts < 0.01) and an explicit 
scheme. One of the main attractions of our scheme is the ability to 
accurately reproduce the forces in poorly resolved flows, i.e. those 
with h > Csts. 

For non-linear drag laws expressions equivalent to equations 
and cannot be expressed in terms of elementary functions, 
except in the case — dig — VP/pg = 0 (although see|Laibe & 


|Price| ( [ 20 T 4 ] ) for the case of quadratic drag laws). In this case the 
only options are explicit integration or an operator split approach. 
In protoplanetary discs the transition to non-linear drag regimes 
occurs for particle sizes a > 6.6 {M/Mmmsn)-\R/ cm, 

where M/Mmmsn is the ratio of the disc mass to the minimum 
mass of the Solar nebula ( [Johansen et al.|20T4] ). Since these parti¬ 
cles have St > 1, the time-steps are limited by the orbital velocity 
and At <C ts, which means that explicit integration is suitable. 

Using analytical solutions means that the time-steps are not 
limited by the stopping time, and we can use the standard limits 
from the gas time-step. We find that a Courant-like condition is 
sufficient to make sure that does not change too much over a 
time-step. We use a signal velocity 


Usig = max(|vd - v^|,Cs), 


(18) 


and set the time-step via dtc = Xi^/'^sig, with xi ~ 0.1. For the 
sound-speed, Cs, we use the sound speed interpolated at the loca¬ 
tion of dust particles. Since the Courant condition is proportional to 
resolution, the dust velocity converges for /i ^ 0, even for ts 0. 
Additionally, we limit the time-step based on the gravitational ac¬ 
celeration, 6ta = X 2 v^/i/|a|. The limits have been sufficient in all 


^ The motivation for this is just that we want to ensure we can evaluate Vg 
at the location of the dust. We show the full problem using interpolation and 
SPH forces in sectionjS^ 
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Figure 1. Convergence of the drift schemes with the number of steps per 
period, N, in the forced oscillation (circles) and sound wave propagation 
(triangles) tests, as measured by the L 2 error norm. The hlled symbols show 
the second order scheme (equationfl^ and the hollow ones show the first or¬ 
der scheme (equation [T^. When the spatial dependence is included (sound 
wave test), the accuracy is slightly degraded, but it is clear that equation [Tb] 
remains second order. 

tests presented, however it may be necessary to include an addi¬ 
tional limiter in the presence of strong shocks (|Durier & Dalla Vec-| 
|chia|2012| l. 


3 NUMERICAL VALIDATION 

In order to validate the scheme, we present several simple 
tests based upon the DUSTYBOX, DUSTYWAVE, settling, and 
DUSTYSHOCK tests described in |Laibe & Price] ( [2011] ) and ( [Loren- 1 
[Aguilar & Bate|20f^ , evaluated in the test particle limit. Using the 
test particle limit means it is possible to compare the simulations 
to semi-analytical solutions, of which few are available for general 
problems. Additionally, we consider multidimensional test prob¬ 
lems, which are more stringent tests of the codes robustness and 
performance in problems that are closer to real physical scenarios. 

3.1 DUSTYBOX for constant drag coefficients 

The DUSTYBOX test is conducted in a 3D periodic box of uniform 
density at a resolution of 20^ gas particles placed on a Cartesian 
grid. We use 20^ dust particles, although since the motion of the 
dust particles are independent a single dust particle could be con¬ 
sidered. The dust particles are placed on a grid offset from the gas. 
The initial dust velocity is taken to be Vcz = 1 in code units, and the 
gas velocity Vg = O.The particles are then evolved using the ana¬ 
lytical solution time-steps for different drag regimes, which is pos¬ 
sible for the non-linear drag laws since there are no external forces. 
In each case the drag coefficients are chosen such that ts = 1 at 
t = 0. 

This test should be trivially passed by our time-stepping 
scheme since the velocity update is exact in this case. Neverthe¬ 
less, this verification is important since some early schemes strug¬ 
gle with even this simple test, although generally this can be solved 
with appropriate kernels ( [Monaghan & Kocharyan|I9^ [Laibe &[ 
[Price[201^ . As can be seen in Fig.|^ the scheme successfully cap¬ 
tures the drag forces in each regime, achieving a force error of less 



Figure 2. Time evolution of the relative velocity between the gas and dust, 
Av = — \-g, for power-law drag forces, oc — | Av|^, with varying 

index k. In each case the coefficients are set such that the stopping time 
ts = latf = 0. The SPH solution is given by the circles, while the 
analytical solution is given by the solid lines. The scheme is able to correctly 
reproduce the force, producing a peak velocity error of 3 x 10“^. 



0.0 0.2 0.4 0.6 0.8 1.0 


X 

Figure 3. Dustywave linear sound-wave test for various stopping times, 
ts. The red lines show the exact solution for the gas (solid) and dust 
(dashed), while the black points and circles show the numerical solution 
for the gas and dust respectively, using 128 particles per phase. The results 
show good agreement between the numerical and exact solution. The L 2 
error norms are the same for all stopping times, but larger than the idealized 
case using the exact gas velocity for the force. The explanation is that the 
SPH force evaluation is the dominant source of error. 

than 0.5 per cent, which verifies the accuracy of the interpolation 
scheme. The results are independent of the ratio of stopping time, 
ts, to gas time-steps, h/cs. Similar results are found using an ex¬ 
plicit time-stepping scheme, albeit at significantly higher computa¬ 
tional cost when Csts h. 

3.2 DUSTYWAVE 

The DUSTYWAVE test ( [Laibe & Price|20lf] ) is essentially the trav¬ 
elling wave test that we used to verify the order of the time-stepping 
scheme in section [2.2[ but using live SPH particles and interpola¬ 
tion to calculate the velocity at the location of the dust particles. 
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We have computed the test in one dimension, using 128 particles 
per phase. We switched off the artificial viscosity since our aim is 
to study the motion of dust particles rather than viscous dissipation, 
which depends sensitively on the artificial viscosity implementa¬ 
tion (see, for example, [Cullen & Dehnen|2010| ). 

In addition to the velocity perturbation 5v/cs — 10“^, we 
set up a density perturbation 5p/p — 10“^ by perturbing the ini¬ 
tial particle positions according to the integrated density along the 
wave. We have used an adiabatic equation of state, with 7 = 5/3, 
and set the entropic function A = to be a constant in order 

to produce Cs = 1 in the unperturbed flow. In each test we have set 
the stopping time ts to be a constant. 

The results of the live test are shown in Fig. which show 
good agreement between the analytical solution and the numer¬ 
ical results. Although not shown, we see similarly good agree¬ 
ment for arbitrarily small stopping times, which shows that our 
numerical scheme successfully captures the evolution of a time- 
dependent fiow, even when h < Cgts and At > which occurs 
for ts < 0.01. In these tests, we find that the L 2 error norm for 
the dust velocity is essentially equal to the error norm for the gas 
and is independent of stopping time. This is due to the dominant 
source of error coming from the SPH force evaluation giving rise 
to a slight bias in the gas velocity and effective sound speed. We 
note that the good agreement in the limit of small stopping times is 
aided by the fact that the DUSTYWAVE test is much less stringent 
in the test particle limit, because unlike the full two-fiuid case the 
sound speed is not modified by the presence of dust, and numerical 
dissipation does not occur, which can be an issue for h < Cstg and 
large dust-to-gas ratios (jLaibe & Price|201 4 [[Loren- Aguilar & Batej 


3.3 Settling in the Epstein regime 


This test involves the settling of particles under the action of gravity 
in a stratified density medium, and is useful for testing the scheme 
in the presence of a varying drag coefficient. For this test we mimic 
a protoplanetary disc in one dimension by applying an external 
force ag — ad — —^^z. We fix the sound speed, Cs, and orbital 
frequency, Q, to 1. The gas is supported against gravity via pressure 
forces, which gives rise to a Gaussian density structure for the gas. 


(19) 


Loren- 


p{z) =poexp , 

where iif = Cs/Q = 1 and po = S/\/27ri/. Following 
[Aguilar & Bate[ { [2014[ ), we use 100 gas particles initially spread 
over the range —2 to 2 , with an initially constant density p(t — 
0) = 1, which implies E = 4, and relax the gas until it reaches 
equilibrium. Once the gas has reached equilibrium we introduced 
100 dust particles, spread uniformly over the range —2 to 2 as be¬ 
fore. We mimic the Epstein regime by setting the stopping time 
ts — ml {Kp), and start the dust particles from rest. 

In Fig. 1^ we show the motion of a particle initially at z = 
—2 for a range of Kjm. In all cases, the results show excellent 
agreement between the numerical solution and the exact solution, 
which for strong coupling tends quickly to the terminal velocity. 
We contrast our results with those reported by [Loren- Aguilar &[ 
[Bate[ ( [2014[ ), who found spurious oscillations in the velocity for 
Kjm ^ 100 at the same resolution. For Kjm — 100, the stop¬ 
ping time is ~ 0.05 at z = ±2 and the smoothing length h ^ 0.2. 
For reasonable values of the Courant parameter, xi = 0.1, we find 
that the time-step. At, is of order the stopping time. We suggest the 
oscillations in their scheme may be due to the incorrect approach 


to the short-friction time limit for At > ts, rather than low number 
of particles. This is consistent with the disappearance of the oscil¬ 
lations at higher resolution, since the Courant condition enforces 
smaller time-steps and thus At < ts. If our explanation is correct, 
then at every resolution there should be minimum stopping time be¬ 
low which their scheme produces oscillations. This test shows the 
benefit of using the full equation of motion for the velocity differ¬ 
ence, rather than an operator splitting scheme. 


3.4 DUSTYSHOCK 


For the DUSTYSHOCK test, we set up the initial conditions in 3D. 
The SPH particles are initially placed on a lattice with a spacing 
of 240 particles per unit length in the high-density region (x < 0), 
152 particles per unit length in the low-density region (x > 0) and 
64 neighbour particles. The masses and internal energies are set to 
give p = 1 and P = 1 for a; < 0 and p — 0.25 and P = 0.1795 
and the adiabatic index 7 = 5 / 8 . 

The dust particles are set up using the same density and num¬ 
ber of particles as the gas. Even at arbitrarily low dust resolution 
the correct velocities are recovered, but an accurate density esti¬ 
mate requires a reasonable number of particles. An equal number 
of gas and dust particles has the advantage that the density in both 
phases experiences the same broadening, making a direct compari¬ 
son simple. Similarly to the DUSTYBOX test, the dust particles are 
placed on a lattice which is offset from the gas. The dust particles 
are then evolved using a linear drag law with a constant stopping 
time, ts. 

In the dust free case, the analytical solution for the gas is well 
known, which equally applies when dust is present in the test par¬ 
ticle limit, (see for example |Rasio & Shapiro|19M) . In addition to 
the analytical solution to the gas, it is possible to compute the solu¬ 
tion for dust via direct integration. For particles that initially have 
X > 0, the solution can be expressed parametrically in terms of the 
initial position of each dust particle, xo = x(t = 0) and the time 
at which the shock reaches the particle to = xq/vs. For linear drag 
forces with constant coefficients this gives 


v{xo,t) 


0, 

Vm(l - exp[-(t - to)/ts]), 


and 


x(xo, t) = Xo + 


{ 


0, 

Vm{t - to) 


v(xo, t)ts 


Xo ^ Vst, 
Xo < Vst, 


( 20 ) 


Xo^Vst, 

XO < Vst, ^ ^ 


where Vs is the shock velocity and Vm is the velocity of the contact 
discontinuity. Since for a stationary shock V • (pv) = 0, the dust 
density can then be written in terms of the pre- and post-shock ve¬ 
locity, p{xo,t) = prVs/{vs — u(xo, t)), whcrc pr is the initial dust 
density for x > 0. For non-linear forces similar expressions can be 
derived straight-forwardly. 

The results of the DUSTYSHOCK test are shown in Fig.|^for 
a linear drag law with stopping times, ts = 0.1 (left) and 0.01 
(right). There is good agreement between the analytical solution 
and the SPH solution. However, for moderate values of the stop¬ 
ping time, ts ~ 0.1, it can be seen that the dust velocity is slightly 
larger in the simulations than the analytical solution. This arises 
due to the broadening of the shock and is most severe for the short¬ 
est stopping times. However, since the particles with the smallest 
stopping times reach the terminal velocity quickly the error in the 
position is bounded. For h/ts < |Av| we find the positional er¬ 
ror is approximately 0.25/i, where h is the gas smoothing length, 
but for longer stopping times the error is smaller. For the smooth 
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Figure 4. Settling test for a range of drag coefficients, showing the evolution of the velocity of a particle starting from rest at 2 : = 2. The solid line shows 
the terminal velocity, while the dashed line show the exact solution for the particle velocity and the circles denote the numerical solution, which uses 100 
particles per phase. We see excellent agreement between the gas and dust and even for K/m > 100 we do not see signs of the spurious oscillations reported 
by [Loren- Aguilar & Bate|{20T^ . 


post-shock flow in the DUSTYBOX problem the post-shock density 
is recovered. However, for more complex flows the positional er¬ 
ror means that SPH simulations of gas and dust are limited to first 
order accuracy. 

3.4.1 Velocity noise and the consequences of too little artificial 
viscosity 

In SPH, shocks compress the particle distribution anisotropically, 
which leads to re-meshing noise as particles re-arrange themselves 
to maintain a uniform distribution ( |Price|[20T^ . Using smoother 
kernels, such as the Wendland or quintic spline, can reduce this 
noise ( |Price|[20T^ jPehnen & Aly|[2012| ). Since the re-meshing 
noise is driven by pressure forces, the dust particles do not undergo 
re-meshing, although the drag forces introduce a velocity disper¬ 
sion into the dust. To measure the velocity dispersion in the shock 
tube tests it is necessary to subtract off the gradient of velocity 
since the broadening of the shock introduces a slight bias into the 
velocity of the dust particles. To this end we use the linear-exact 
SPH gradient estimator ( [Price |20T^ , which results in an estimate 
of cr^ = 3 X 10“^ Cs in the post-shock gas. For the dust particles 


with ts — 0.01 the velocity dispersion is the same as the gas. For 
more weakly coupled dust the velocity dispersion is lower, with 

Gv < 10“^Cs forts = 0.1. 

While these noise levels are considerably lower than the veloc¬ 
ity dispersion typically measured in simulations of self-gravitating 
protoplanetary discs, we do And that the noise is sensitive to the 
value of the artificial viscosity parameter, a. As described above, 
early simulations of self-gravitating discs, including the first with 
dust, employed a = 0.1 instead of the usually employed a = 1. 
As a result, the simulations fail to generate sufficient entropy in 
shocks which gives rise to oscillations in the velocity and density 
behind the shock. We illustrate this in Fig.|^ which also shows the 
considerable noise in the dust velocity driven by these oscillations. 
We note that in this case, subtracting off the gradient produces a 
considerably lower estimate for the velocity dispersion in the gas 
(gv ~ 10“^Cs) than if we use the difference between the simula¬ 
tion and analytical solution for the estimation (gv ~ O.lcs). For 
the dust the difference between the estimators is smaller, since the 
oscillatory structure is less well defined in the dust. 

Since the shocks in protoplanetary discs are weak, we have 
measured the velocity dispersion generated by the re-meshing noise 
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Figure 5. Shock-tube tests for combined gas (black) and dust (red). The tests are calculated in 3D with stopping times ts =0.1 (left) and ts = 10“^ (right). 
The analytical solution for the gas (blue) lines and the dust initially to right of the contact discontinuity (green) are also shown. The choice of a cubic spline 
kernel and 64 neighbours gives rise to significant re-meshing noise in the gas velocity. However, the interpolated gas velocity at the location of the dust is 
considerably smoother, resulting in significantly less noise in the dust particles. Typically the velocity dispersion in the dust is less than 1 per cent of the sound 
speed, compared with 5 to 10 per cent for the gas. 
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Figure 6. Velocity in the shock-tube test with ts =0.1 (left) and ts = 0.01 (right) when a low viscosity parameter a = 0.1 is used. The phenomenon of 
post-shock ringing in the gas is well known. While this translates into post-shock noise in the dust velocity, the amplitude is much smaller as the interpolation 
averages over the gas velocity. 
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Figure 7. Velocity dispersion in the post-shock gas in shock tube tests with 
varying Mach number, Ai. The different symbols (triangles, squares, pen¬ 
tagons) refer to dust particles with differing stopping times, and the corre¬ 
sponding velocity dispersion in the gas. The hlled symbols refer to standard 
simulations, whereas the hollow symbols show simulations where a low ar¬ 
tificial viscosity, a = 0.12, was used. The points show the median velocity 
dispersion measured and the error bars show the 25th and 75th percentiles. 
We show Ai — Ion the x-axis to highlight the low velocity noise for weak 
shocks, as expected in protoplanetary discs. 


for Mach numbers between 1.03 and 1.44, which are shown in Fig. 
|7] In this range, we find Gv oc Ad — 1, reflecting the relationship 
between the compression in the post-shock gas and the velocity 
noise in the gas. The low velocity dispersion introduced in the dust 
particles by weak shocks Gv < 10 shows that even though 
a = 0.1 under-produces entropy in shocks, planar shocks are un¬ 
likely to be the dominant factor determining the velocity dispersion 
of dust in simulations of protoplanetary discs. In the following sec¬ 
tions, we consider multidimensional tests to explore whether the 
level of noise produced in the shock tube tests is typical of more 
general problems. 

We note that the level of velocity noise reported is implemen¬ 
tation dependent and refiects the levels of noise expected in a com¬ 
monly used SPH implementation. In particular, the results will be 
sensitive to the choice of kernel, since smoother kernels that use a 
larger number of neighbours, such as the Quintic Spline or Wend- 
land kernel, produce lower levels of noise in the gas velocity 
( |Price|20f^|Dehnen & Aly|2012| l. In cases where numerical noise 
from shocks needs to be minimised, these kernels are highly recom¬ 
mended since the dust particles lack the ability to regularize their 
particle distribution (see also section [33] ). 


3.5 Shocks in 2D 


The high symmetry of the Sod shock tube problem means that the 
uniformity of the initial conditions is preserved. However, in real 
applications the high symmetry of the initial conditions is broken 
as the system evolves. In order to show the effect that this has on 
the dust distribution we have conducted the implosion test ( |Hui,| 
[Li & Li|199^, which we have modified to use periodic boundary 
conditions ( ^jacki et al.|2012| ). For the initial conditions we con¬ 
sider a two-dimensional box with —0.15 < x^y < 0.15. The gas 
is initially at rest with density and pressure are p = 1 and P = 1 , 
except for X -\- y < —0.15, where p — 0.125 and P = 0.14, and 
the adiabatic index 7 = 1.4. We have set up the test using equal 


mass particles and Nx — Ny — 600 particles in the high density 
region. The particles have been on a hexagonal lattice and we use 
14 neighbour particles in 2D (equivalent to 39.4 neighbours in 3D). 
We add dust to this problem assuming a uniform initial density, us¬ 
ing Nxg = Nyg = 600 particles placed on a hexagonal lattice 
offset from the lattice of gas particles. The stopping time is taken 
to be constant, ts — 0 . 1 . 

The shocks launched by the initial discontinuities interact in 
the low density region and should drive a jet of gas along the sym¬ 
metry axis X = y. The interaction of this jet with multiple shocks 
leads to Richtmyer-Meshkov instabilities. Vanilla SPH, as imple¬ 
mented in GAD GET- 2, is known to struggle with such instabilities 
( jAgertz et al.|2007| ). Although we are more interested in how the 
dust evolves in a problem without a high degree of symmetry, we 
have included an artificial conductivity as described by |Price|(2012| ) 
and a viscosity limiter ( [Cullen & Dehnen|2010| l to help reduce the 
damping of instabilities. 

We have also run the test using the grid-code FARGO ( [Masset 
2000 1 , which is based upon the ZEUS algorithm ( [Stone & Norman 
1992 1 to provide a reference solution for the dust. Dust has been 
included in FARGO using a finite volume approach in a similar 
way to jZhu et “^ ( [20121 ), using explicit time-stepping for the dust 
acceleration. Strictly a finite volume approach is not valid for this 
problem, since the dust velocity becomes multivalued in the region 
x,y < 0. While this results in density artefacts that arise on the grid 
scale, this does not degrade the solution elsewhere, which can be 
confirmed by comparison with the SPH simulations. The FARGO 
simulations were run using 600 x 600 cells. 

In Fig. [^ we show the density at t = 0.7, after approximately 
a sound-crossing time. In the simulations in the top row, the large 
scale density structure is reasonably well reproduced, and the com¬ 
bination of a viscosity limiter with artificial conduction allows the 
growth of instabilities. However, the instabilities remain weak due 
to large amounts of noise in the density and velocity on scales of a 
smoothing length. The noise is particularly evident in the vorticity, 
which is largely featureless. The failing of standard SPH to repro¬ 
duce the vorticity was noted by jSijacki et al.| ( [2012| ), and we find 
that while using a viscosity limiter and artificial conductivity helps 
the instabilities grow, it also results in a more noisy solution. The 
larger noise in the vorticity is associated with the lower viscosity - 
using a fixed a = 1 not only damps the instabilities but also damps 
the velocity noise and helps maintain a regular particle distribution. 
The noisy gas velocity translates into considerable noise in the dust 
density on scales close to a smoothing length, even in regions that 
are not affected by the instabilities. The difference in the noise in 
the dust and gas densities provides a clear illustration of the impor¬ 
tance of the self-regulation in SPH in producing an accurate solu¬ 
tion, since the dust has no mechanism by which it can reduce the 
noise in the density distribution once it has been introduced. 

Since the higher density noise in the dust is associated with 
the fact that the dust particle distribution is less regular, the noise 
could be decreased by increasing the number of neighbours. For 
a standard cubic spline kernel particle, the pairing instability pre¬ 
vents the use of many more neighbours in the gas ( [Price|[20T2| ). 
Since the drag forces have been included using averages of the 
gas velocity the dust particle should not be affected by the pair¬ 
ing instability. In this case it makes more sense to spend additional 
computational time to ensure the gas dynamics is computed accu¬ 
rately. To demonstrate this, we recompute the test using a modem 
SPH implementation that uses a Wendland kernel, 50 neighbours, 
and integral based derivatives. The implementation is essentially 
identical to the formulation of |Rossw^ ( |20I4j ). Increasing the 
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Figure 8. Gas and dust density (left & middle) and gas vorticity (right) in the 2D shock problem at f = 0.7. The top row shows the results for SPH using 
artificial conductivity and the [Cullen & Dehn'^ ( |2010) viscosity switch along with the cubic spline kernel. In the middle row a Wendland kernel and 50 
neighbours were used, along with integral-based derivatives. In the bottom row, for comparison we show the results of the test problem as calculated by the 
grid-code FARGO, with the dust treated using a finite-volume approach. The noise in the gas density and velocity in the top panels has prevented the growth 
of Richtmyer-Meshkov instabilities, and is responsible for a large amount of noise in the dust density. 


number of neighbours results in an additional computational cost of 
a factor of three. Since the integral based derivatives are required 
for the [Cullen & Dehn^ ( |20I0| ) viscosity switch, they can also be 
used for the force calculation at negligible extra cost, although we 
note that this represents a 20 per cent increase over simulations in 
which the viscosity parameter is kept at a fixed value of a = 1. 

The results using this scheme are shown in the middle row of 
Fig|^ The effect of these modernizations on the gas density is rel¬ 


atively minor, resulting in sharper shocks and stronger Richtmyer- 
Meshkov plumes. The sharper shocks are a result of the reduced 
noise and give the appearance of higher resolution, although the 
same resolution has been used, albeit at considerable extra compu¬ 
tational cost. However, the improvement that this makes on the dust 
density is staggering. The density noise in smooth parts of the fiow 
is significantly reduced and vortices in the gas appear that are de¬ 
void of dust. These results are in good agreement with the FARGO 
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calculations, shown in the bottom panel of Fig. which exhibit a 
higher effective resolution. The main difference arise in the region 
x,y < 0, reflecting the both the higher resolution of the FARGO 
simulations in this region, along with the issues related to the dust 
velocity being multi-valued in this region. 

The vast improvement in the dust density distribution with the 
state of the art SPH calculation reflects the fact that reasonable 
agreement in the gas density can be achieved even when the gas 
velocity structure is dominated by noise. Since the gas velocity de¬ 
termines the dust dynamics, the dust density reflects the accuracy 
achieved in the velocity structure. The modem SPH implementa¬ 
tion clearly does a much better job at reproducing the gas velocity, 
which also results in a much better reproduction of the instabili¬ 
ties. As this test shows, an accurate velocity structure is essential 
for reproducing the dust dynamics and convergence in the velocity 
structure can be considerably more difficult to attain. 

3.6 Static disc problem 

In addition to the simple veriflcation tests above, we present a test 
for the code in a mode of operation that is close to the conditions 
present in protoplanetary discs. For this, we consider the motion of 
an isothermal gas in a flxed potential, consisting of an axisymmet- 
ric background potential with an imposed spiral perturbation. For a 
logarithmic potential the solution can be calculated along stream¬ 
lines that orbit at a given average radius ( [Roberts |1 969 1 [Shu, Mil-~| 
jione & Roberts] 197^. In partic ular, we follow the method of solu- 
tion given in jGittins & Clark^ ( |2004| ). We use a logarithmic spiral 
perturbing potential, 14, with a constant pitch angle, i, given by 

14 = AoRex.p{-esR) cos(x), (22) 

where 

777 

X =-: In(esi^) — m(0 — ^pt). (23) 

tanz 

For the number of spiral arms, m, we use m — 2, the pitch angle, 
sinz = 0.1 and the pattern speed at which the spiral potential ro¬ 
tates, Qp = 0.522. The constants Aq and Cg set the strength of the 
potential and the length over which it decays, and R is the cylindri¬ 
cal radius. The spiral perturbation strength can be set relative to the 
background potential, we take Cg = 0.85 and F = 0.05 at = 1, 
where 

FIR) = |VK|/|VM,| = (24) 

For the background potential, a Keplerian velocity Q(i?) = 
^GM/R^ could be used. However, the streamline solutions can¬ 
not be found in regions close to corotation, where the background 
velocity is subsonic, R(Q — Qp) sin Since the spiral modes 

present in protoplanetary discs are close to corotation it is not pos¬ 
sible to directly probe the regime present in protoplanetary discs 
using the streamline solutions ( jCossins, Lodato & Clarke|200^ . In¬ 
stead, we use a flat rotation curve in which the streamline solutions 
have been well studied in relation to spiral galaxies. Nevertheless, 
this allows us to test how the code handles shocks in a shearing 
environment. For the background potential we use the form and 
parameters from jGittins & Clarkej ( [2004] ), which we scale to di¬ 
mensionless values. The velocity is given by 

v(R) = FbebRe^p{-ebR) + 1 - exp(-edi^), (25) 

where we set the maximum rotation velocity Umax from v{R) = 1 
at = l.The parameters Fb, and deflne the strength and size of 
the galactic bulge, which we neglect by taking i4 = 0. The scale 


of the disc is set via = 17/3 and the sound speed of the gas is 
= 27.5. 

The streamline solutions are calculated in coordinates that de¬ 
note the position along, and between, 77, the spiral arms. Choos¬ 
ing ?7 = TT — X means that 77 = 2mT for integer n at the potential 
minima. Choosing ^ to be perpendicular everywhere to 77 gives 

Tfl 

rj = -In(esi^) + m{0 — Qpt) + tt (26) 

tanz 

777 

^ = -m\n(egR) +- AO - Qpt). (27) 

tanz 

The streamline solution can be calculated at a given radius un¬ 
der the approximation that the local angular velocity, epicyclic fre¬ 
quency, perturbation strength and sound speed are constant along 
the streamline. The streamlines are found by choosing a value r]g 
for the sonic point and integrating forward and backwards along 
the streamline, whilst looking for a shock. Once a shock has been 
found 77s is varied until a period solution is found. For more details 
on the exact method, see Appendix A of jGittins & Cl^ ( |2004l ). 
In addition to the equations for the velocity along and perpendic¬ 
ular to the streamlines, Vp and u^, we solve for the position about 
the disc by integrating 


along a streamline ( [Roberts jl 969 j ). 

Once the streamline solutions have been calculated, the full 
structure of the gas can be built by interpolating between them. 
However, the periodicity of the streamlines is only approximate, 
with the initial and flnal radii differing by roughly 1 per cent after 
the gas has travelled through an angle of 27r jm. For this reason 
instead of computing streamlines that give velocities that are peri¬ 
odic over A77 = 27r, we look for solutions that are periodic over 
AO = 27r/77z. This is motivated by the need for streamlines that do 
not overlap themselves in the interpolation process. 

We now simulate the evolution of dust particles that evolve 
under the action of the spiral gravitational potential as well as drag 
forces exerted by gas whose density and kinematics is given by the 
analytic gas streamline solutions. Henceforth the dust properties 
derived from these simulations are described as ‘dust streamline 
solutions’ and these are to be compared with the properties of the 
dust population obtained when the dust is instead introduced into a 
‘live’ SPH calculation and the dust-gas drag computed as described 
above. For the ‘dust streamline solutions’, up to 10^ particles are 
placed on a grid with an initial velocity taken to be the unperturbed 
background velocity. The particles are then evolved using a fourth 
order Runge-Kutta-Fehlberg integrator. Once all the dust positions 
are known, the dust density in any snapshot is then calculated using 
an SPH kernel sum. 

The SPH simulations are conducted both in 2D and 3D. In 
the 3D simulations, periodic boundary conditions are used in the 
^-direction to allow direct comparison between the 2D and 3D 
results. This avoids the additional complications associated with 
stratifled discs and the settling of dust into the mid-plane which 
arises within them. While the difference between velocity distri¬ 
butions in 2D and 3D is an interesting issue in its own right, it 
complicates the comparison to the analytical solution. 

The SPH particles are set up using a uniform density glass, 
and the gas is evolved under the static potential until they reach a 
steady-state, roughly 35 orbits at = 1. Once the gas has reached 
a steady-state the dust particles are introduced and the system is 
evolved for a further 35 orbits. Tests conducted in which the per¬ 
turbing potential was slowly applied to the gas show no signiflcant 
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Figure 9. Gas (left) and dust (right) density in the isothermal disc for the 3D simulation with h/H = 0.16 at = 1 and Stokes number, St = 1. The dust 
density has been computed using twice the smoothing length of gas, to reduce the noise in the density estimate. The white line shows a gas streamline for 
comparison. The flow is anticlockwise and the density maxima in the dust appears behind the density maxima in the gas. 




Figure 10. Convergence of gas density and velocity perpendicular {vr^) and parallel {v^) to the spiral equipotential lines at = 1. The potential minima 
and maxima are denoted by 77 = 0 and 77 = ± 7 r. The SPH simulations include two 3D simulations with resolution h/H = 0.35, and 0.16 (magenta and 
green lines) and a 2D simulation with h/H = 0.05 (yellow). The black line gives the streamline solution. The SPH simulations converge well to the expected 
density and perpendicular velocity. However, the parallel velocity converges to a velocity that is 1 per cent lower than streamline solution. 


differences to tests in which the potential was switched on imme¬ 
diately. In all results presented the latter approach was taken. The 
uniform density means the inner regions of the disc are less re¬ 
solved than the outer regions; therefore, it is instructive to quote the 
resolution in terms of h/H, where h is the smoothing length and 


H = Cs/Q is the pressure scale-height. For comparison at = 1, 
3 . h/H of 0.35 and 0.17 translates into 2.5 million particles with 
a box-height of Z = 3 and 6 million particles with a box-height 
of 0.75. lnlT> h/H — 0.05 can be achieved with 4 million par¬ 
ticles. The number of particles used here is larger than the ~ 10® 
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Figure 11. Dust properties along a streamline for stokes number St = 1. The black lines denote the gas and dust properties calculated from the streamline 
solution and the results from the SPH simulations are given by the coloured lines. The dust density has been calculated using the equivalent of 50 (bottom left) 
and 400 neighbours (top left) in 3D. 


particles typically used in protoplanetary disc simulations (c.f.|Rice| 
|et al.|2004[|Cossins, Lodato & Clarke|2009[|Meru & Bate|2()12| ), 

which reflects the fact the model disc is thinner H/R r\j 10 ^ than 
protoplanetary discs {H/R ~ 0.1), resulting in the need for con¬ 
siderably more particles to resolve the pressure scale height. 

The density structure in the disc is shown in Fig.[^ for the 3D 
simulation with h/H = 0.16 and St = 1, upon which a typical 
streamline is marked for comparison. Each streamline crosses both 
of the spiral arms once. In calculating the dust density 400 neigh¬ 
bours have been used. The large number of neighbours needed in 
the dust density calculation reflects the fact that no self-regularizing 
forces apply to the dust. For a comparison of how noisy the dust 
density is with the normal ~ 50 neighbours, see Fig. m The flow 
is anti-clockwise and the dust density peaks signiflcantly down¬ 
stream of the shock, since the dust particles respond to the change 
in velocity over a stopping time. 

Fig. shows the gas density and velocity along a stream¬ 
line at = 1, showing that standard SPH is able to adequately 
reproduce the streamline solutions. While the SPH simulations ac¬ 
curately reproduce there is minor discrepancy in since the 
streamlines converge to a velocity that is systematically 1 per cent 
lower than the streamline solution. While pressure support can re¬ 
duce the velocity below its Keplerian value, the disc is isothermal 
and initially has no density gradient. Once the spiral perturbation is 
applied a density gradient does arise, giving an order unity change 
in the density over the size of the disc. However, since the disc 
is cold and the ratio H/R ~ 10“^, pressure support should only 
contribute at the 10“^ level. 

In some cases, artiflcial viscosity can be responsible for pro¬ 
ducing an incorrect disc structure in SPH simulations, via viscous 


angular momentum transport. For example, it has been shown that 
when a constant a = 1 is used in the differentially rotating Gresho- 
Chan vortex problem that simulations do not converge to the cor¬ 
rect inviscid solution, but appear to converge to a different solution 
( |Springel|20I0| ). However, it is possible to show that artiflcial vis¬ 
cosity is not responsible for the velocity difference observed in our 
simulations. Firstly, in our test problem the spiral perturbing po¬ 
tential is responsible for driving the density gradient. Furthermore, 
the tests with low viscosity also show the velocity difference, as do 
tests in which a viscosity switch similar to |Cullen & DehnelilpOlOl ) 
was used. Since the streamline solutions are only approximate, and 
the streamlines only close to within 1 per cent, the streamline solu¬ 
tions themselves may be the origin of the error. While the discrep¬ 
ancy is clearly greater than the noise in the simulations, the agree¬ 
ment with the streamlines is good enough that we can compare the 
dynamics of dust in the simulations to the streamline solutions. 

In Fig. m the properties of dust for a linear drag law with 
St = 1 are shown along the same streamline at = 1. The density 
has been calculated using both 50 and 400 neighbours, showing 
that a large number of neighbours is needed to achieve an accurate 
density estimate. Increasing the number of neighbours has a much 
more dramatic effect on the density estimate than increasing the 
resolution, which reflects that the local density estimate is subject 
to Monte-Carlo noise once the smooth density held present in the 
initial conditions has been sheared away. Conversely, the noise in 
the dust velocity is considerably lower, and can be clearly seen to be 
converging to the streamline solutions as the resolution increases. 

Since the velocity distribution of dust particles is fundamental 
to understanding the growth and evolution of dust grains in proto¬ 
planetary discs, it is important to check that the velocity distribu- 
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Figure 12. Velocity dispersion as a function of Stokes Number for dust par¬ 
ticles in SPH simulations (black) and streamline solution calculations (ma¬ 
genta). The different symbols mark the resolution of the simulations h/H. 
The sampling resolution of the streamline solution is close to the highest 
resolution SPH simulation. For the SPH simulations, hlled and empty sym¬ 
bols mark simulations with standard viscosity and low viscosity respec¬ 
tively. The points mark the median velocity dispersion and the error bars 
denote the 25th and 75th percentile. The simulations at h/H = 0.05 were 
conducted in 2D, with both the higher resolution and lower number of de¬ 
grees of freedom contributing to the lower noise. 


persion in the simulation is dominated by noise since the velocity 
dispersion of the streamline solutions is considerably lower. How¬ 
ever, the velocity dispersion in this problem is orders of magnitude 
lower than that seen in simulations of the time-dependent structures 
in self-grav itating protostellar discs, which have ay ~ Cs ( |Rice| 
|et al.|2004] ). Since the numerical noise is at the 10level in our 
/i/iT = 0.16 3D calculation, the velocity dispersion found in sim¬ 
ulations of dust evolution is more likely to be associated with the 
physics of the problem, (i.e. the fluctuating spiral structure) than 
being a simple artefact of the numerical implementation. 

We have tested how using a modern SPH implementation af¬ 
fects the gas and dust dynamics in the spiral disc test. Although 
these methods help enormously with the two-dimensional shock 
tube test we found only a minor improvement when using integral- 
based gradients in the spiral disc problem. This can be understood 
since in the simple test cases the Wendland kernel and integral gra¬ 
dients produce a very regular particle distribution. However, in the 
static disc problem the shear prevents this from happening, which 
limits the benefits these improvements can provide. Additionally, 
we have run simulations using a |Cullen & Dehnen] p010| ) limiter 
with Qmax = 1-2, and amin = 0. With these parameters we found 
that Gy is 50 per cent higher than when using a fixed a — 1.2, since 
viscosity acts to damp the noise in the particles. However, using a 
low fixed a = 0.1 results in a Gy is several times higher than is 
produced when using the limiter. Clearly using a |Cullen & Dehnen| 
( |2010| ) limiter is preferred when in problems where it is necessary 
to reduce viscosity away from shocks. 


tion is not dominated by numerics. Since in a shearing disc there 
are variations in the velocity at the order of the sound speed, these 
need to be subtracted off in order to measure the velocity disper¬ 
sion. To do this, we fit a third order polynomial to velocities at each 
point in space. We fit the velocity using a smoothing length that 
gives 400 neighbours to reduce the noise. We then subtract off the 
fitted velocity held and calculate the velocity dispersion from the 
residuals, again weighted by the kernel function. The method we 
use to fit the polynomial is identical to that described in |Read &| 
|Hayfleld| ( |20T^ , except we fit a 3rd order polynomial. For more 
details, see Appendix C of [Read & Hayfleld| ( |2012| ). 

The velocity dispersion for the dust particles for a range of 
stopping times and resolutions is shown in Fig.[^ The use of an¬ 
alytical solutions in the time-step means that we can investigate a 
range of Stokes numbers, St < 0.1, that would be inaccessible us¬ 
ing explicit time integration. Since the resolution varies as a func¬ 
tion of radius, we calculate the velocity dispersion at points along a 
streamline diiR— 1, for both the simulations and the streamline so¬ 
lutions. The simulations And normalized velocity dispersions in the 
solid component of as little as Gy/cg — 10“^ at the highest resolu¬ 
tion and St = 1, which is 3 x 10“^ of the orbital speed. However, 
since this level of velocity is larger than that measured from the 
dust particles in the streamline solutions at the same resolution, it 
can still be considered numerical noise. For St = 10“^ even the 
velocity dispersion from the streamline solutions is dominated by 
numerical noise. 

For the worst case scenario, the lowest resolution simulation 
with a = 0.12, the median velocity dispersion in the gas parti¬ 
cles reaches O.lcg. Typically the low viscosity simulations produce 
5 times more velocity noise than those run with a — 1.2. Sim¬ 
ilarly to the shock tube test, the velocity dispersion in the dust 
is typically a factor of 10 smaller than that in the gas. Compari¬ 
son of the velocity dispersion measured from the simulations and 
streamline solutions near h/H — 0.1 shows that the velocity dis- 


4 DISCUSSION 


For a over a decade SPH has been relied on to provide insights 
into the dynamics of gas and dust in self-gravitating discs. Indeed 
to date, most of the global 3D simulations of self-gravitating discs 
and been performed with SPH (although [Holey & Durisen| { |2010| ) 
investigated enrichment of solids in gas giants that are formed by 
gravitational instability using a 3D Eulerian code). It is therefore 
necessary to take a critical look at the fldelity of the method. While 
many code comparisons exist in the case of purely gaseous discs 
( jMeru & Bate|20lf]|2012||Paardekooper, Baruteau & Meruj20II| l, 

the dust implementation has been hitherto untested in this context 
(although, for studies of dust in general see Laibe & Price||20II| 


[M^lAyliffe et al.|2012|[Loren- Aguilar & Bate|2014t.~ 

A key question in the case of self-gravitating discs concerns 
the velocity dispersion of solid material in the disc, since this con¬ 
trols the prospects for gravitational collapse in the dust phase as 
well as for particle growth or destruction. Previous simulations 
( [Rice et aL||2004| |2006| [Gibbons, Rice & Mamatsashvilij|2012'] ) 

found a large velocity dispersion in the dust (of order a sound 
speed), which is unfavourable to either grain growth or gravita¬ 
tional collapse. It is therefore important to test whether this result 
is likely to be a numerical artefact. 

In order to test the reliability of such codes, we have imple¬ 
mented dust particles into the SPH code GAD GET- 2. We have used 
a two fluid approach and limited the study in the test particle limit 
where the feedback from the dust onto the gas can be neglected, 
allowing a direct comparison between simulations and analytical 
models. The code has been designed to accelerate the time-step us¬ 
ing analytical solutions for the drag forces, in order to cope with 
particles that have Stokes number St <C 1. The chosen method 
has the advantage that it captures the full phase-space information 
of the dust, which is essential for evaluating the growth and sur- 
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vival of dust grains in protoplanetary discs. The code can handle 
the full range of Stokes numbers and produces the correct termi¬ 
nal velocity for At ^ tg. We have focussed on extensively testing 
the code and its ability to tackle the questions in the dynamics of 
gas and dust in self-gravitating discs. In many of the tests presented 
we have focussed on the dynamically interesting cases in which the 
stopping time is comparable, or slightly smaller than the dynamical 
time scale. The result is that our time-stepping scheme only really 
saves computational time for the shortest stopping times. However, 
demonstrating that the time-stepping scheme can handle both situ¬ 
ations readily is important for realistic physical problems, such as 
protoplanetary discs in which both short and long stopping times 
are present. 

The code performs well in the standard DUSTYBox and 
DUSTYSHOCK test problems, in which the high symmetry means 
that the orderly distribution of the particles is maintained by the 
symmetry in the problems. We also find good agreement in DUSTY- 
WAVE test, and have shown that our method is able to correctly 
reproduce the velocity of particles in a settling test, even at low 
resolutions and in strong drag regimes. We also run the two dimen¬ 
sional shock-tube or implosion problem ( |Hui, Li & Li|1999} , which 
involves interacting shock waves and is a much more severe test of 
a code’s robustness. We have shown that correct modelling of the 
dust requires accurate modelling of the gas dynamics; Fig. [^con¬ 
tains a striking demonstration that a correct modelling of the gas 
density distribution does not guarantee that the gas dynamics (and 
the density distribution in the dust) is well represented since the 
dust is quite sensitive to noise in the gas velocity. In agreement with 
the results of |Dehnen & Aly| ( |2012| ) in the case of the Gresho-Chan 
vortex test, we find that using a Wendland kernel (which allows 
the number of neighbours to be increased beyond the traditional 
limit imposed by the pairing instability) makes a huge difference 
to the accuracy of the velocity field, as does using integral based 
derivatives. Since grid codes compare much more favourably in the 
two-dimensional shock tube problem, grid-based implementations 
of gas-dust mixtures will not show the same poor performance that 
standard SPH does in this problem. However, it would be inter¬ 
esting to see how one-fiuid SPH codes perform in the presence of 
interacting shocks. 

We also tested the code on a problem that includes spiral 
shocks in a shearing environment, which is applicable to shocks in 
self-gravitating discs, or structures that form from planet-disc in¬ 
teractions. We find that the code is able to reproduce accurately the 
motion of the dust particles, giving the expected solutions. How¬ 
ever, in general this is less true for dust density. While the gas 
density in SPH calculations is accurate to the 1 per cent level for 
roughly 50 neighbours, we have found that even with 400 neigh¬ 
bours 10 per cent accuracy in the dust density is not achieved. The 
difference arises because the SPH particle distribution is continu¬ 
ally regularized by pressure forces, but the dust distribution is not. 
Once the disc has evolved for a sufficient length of time that the ini¬ 
tial structure of the disc has been sheared away, the dust particles 
then randomly sample the density distribution. This means the den¬ 
sity error, Up oc , for the dust particles while Up oc 

for the gas. Our results agree with the findings of|Zhu, Hernquist &| 
|Li|( |20T4l ), who find that roughly 1000 neighbours are required to 
calculate the density with an accuracy of 10 per cent when particles 
are placed randomly in a volume. The reason this density noise is 
not seen in the shock tube problems is that the turbulence has not 
sufficiently deteriorated the initial particle distribution. 

Since the dust density does not enter the equations of motion 
within the test particle limit, the error does not affect the ability of 


the code to calculate the motion of dust particles. The low velocity 
dispersion of the dust in the steady disc problem ctv < 10ver¬ 
ifies this. The good performance of the code in the static disc test 
problem is promising, producing much lower velocity dispersion 
than seen hitherto in simulations of self-gravitating discs. However, 
there are some caveats. Firstly, using a low viscosity parameter, 
a — 0.1, fails to generate enough thermal energy at shocks, which 
results in post-shock oscillations in the gas when an adiabatic equa¬ 
tion of state is used. As the static disc test case is isothermal the 
energy generation requirements do not apply, which may help to 
reduce the noise generated. However, in all our test cases the am¬ 
plitude of the velocity dispersion is smaller than the velocity jump 
across the shock Uv/cs < M, setting a natural limit on the amount 
of noise introduced by shocks. It remains to be seen how the dusts 
dynamical state will be affected by these various code choices in 
the self-gravitating case and we leave this exploration (which will 
help to distinguish between physical and numerical effects) to a 
future paper. 

Furthermore, when the dust is sufficiently dense to affect the 
dynamics of the gas, reducing the density error, which can be of 
order unity, will be important for producing reliable results. This is 
one area in which using a single-fiuid approach would help, since 
the dust properties are evaluated at the location of the SPH parti- 
cles and the ap^ oc noise does not apply. However, in the 

presence of velocity noise the solution obtained by one-fluid ap¬ 
proaches may also be deteriorated. Furthermore, in problems where 
the full phase space information is required, using a two-fiuid ap¬ 
proach with a large number of dust particle neighbours may be nec¬ 
essary to ensure an accurate force estimate. Indeed, similar results 
were reported by |Laibe & Price|f2012| ), who find significant noise 
in the dust density in a protoplanetary disc simulation, which led to 
artificial clumping. They found that using the gas smoothing length 
partially resolves the issue, since the gas has lower density and 
a larger smoothing length. We suspect that both error in the dust 
density and requirement that the gas is resolved on length scales 
present in the dust distribution are important for resolving the issue 
( [Ayliffe et al.|2012| ). 

Reducing the noise in the dust density distribution is even 
more important when dust self-gravity is important, since spurious 
over-dense clumps may collapse and feed back on the gas phase 
fragmentation. When self-gravity is important it is likely to be es¬ 
sential that a large softening length is used. Therefore, care must 
be taken in simulations of fragmentation in the dust layers of self- 
gravitating discs to avoid numerical artefacts. This is likely to ap¬ 
ply to all particle based implementations of dust, independent of 
whether the gas is calculated using SPH or with a grid code. 


5 CONCLUSIONS 

We have implemented dust particles in the Smoothed Particle Hy¬ 
drodynamics code gadget-2. Similarly to [Loren- Aguilar & Bate] 
( |2014j ) we avoid the need for small time-steps when the dust and 
gas are tightly coupled. However, unlike previous methods, our 
time stepping scheme produces the correct terminal velocity even 
when the time-step is much greater than the stopping time. We have 
applied a number of simple test problems that verify the code pro¬ 
duces the expected results. The tests show the importance of pro¬ 
ducing an accurate velocity structure for capturing the dust dynam¬ 
ics, since a noisy velocity field introduces noise into the dust den¬ 
sity that persists for the duration of the simulation. 

We have used a rigidly rotating potential to test the dynam- 
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ics of dust in spiral discs, and find that as long as the gas is 
well resolved numerical noise contributes a velocity dispersion 
Gv ~ 10“^ to 10“^Cs, below the levels found in simulations of 
self-gravitating protoplanetary discs. We find the dust density in 
disc problems is more difficult to estimate as shear destroys the 
high accuracy present in the initial density field. This means that 
the dust density is subject to Monte-Carlo noise and more than 400 
neighbours are required for an accurate density estimate. While er¬ 
rors in the dust density do not affect the tests presented here, which 
are calculated in the test particle limit, when the dust density is high 
enough to be dynamically important, either through drag forces or 
self-gravity, care must be taken to ensure accurate results. 

From these numerical experiments we have shown that for 
simulations in which the length scales in the gas are well resolved, 
h/H ~ 0.1, the dynamics of dust can be accurately modelled. 
These tests have been conducted in the test particle limit and in 
relation to SPH, but we expect them to hold more generally with 
the caveat that if the dust reaches high enough densities to affect 
the gas dynamics then a large number of neighbours is needed to 
ensure robust results. As long as care is taken to ensure the density 
noise is controlled, then simulations should be able to tackle some 
of the key questions related to the physics of dusty gasses. 
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APPENDIX A: KERNELS & INTERPOLATION 

The evaluation of drag forces in SPH requires interpolation of the gas prop¬ 
erties to the location of the dust particles and vice versa when the feedback 
on the dust is important. The most straight-forward way that this can be 
done is via SPH interpolation at the location of the dust particles, giving a 
force 

Fd = -V—i^id(vd-Vi)W(rid,ftd), (Al) 

i 

which is the approach taken by [Rice et al.H2004} and used throughout this 
paper. However, if this approach is used for the velocity difference then the 
expression does not conserve angular momentum explicitly jMonaghan &| 
|Kocharyan| 1995 [[Laibe & Price|2012^ . To solve this they propose projecting 
the force been pairs of particles along the line between them, which gives 

Ed = (A2) 

where a is the number of dimensions and = 'Vd — Equation |A2[ 
approximates the force as the average of a series of forces acting in dif¬ 
ferent directions. This means that even if the velocity field is uniform, the 
distribution of particles affects the force. Additionally, even if the particle 
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Figure Al. Accuracy of different formulations of the drag law for the 
DUSTYBOX test, for the formulations given by equations [A^ and [A^ Both 
formulations produce similar results as long as appropriate kernels are used. 
While increasing the number of particles does not improve the accuracy of 
the drag forces, increasing the number of neighbours does, since this im¬ 
proves the local estimate of the gas properties. 


distribution is regular a large number of neighbours is needed to evaluate 
the force accurately. |Laibe & Price] suggest resolving this problem 

by using a double hump kernel that has the form W(q) oc q‘^G{q), where 
G(q) approximates a Gaussian. This essentially removes the contribution 
of particles closest to the dust particle. Since it is these particles that drive 
the bias when a centrally peaked kernel is used, down weighting them is 
effective in removing the bias. For the original form (equationthe bias 
is not present and a kernel that approximates a Gaussian should be used. 

We illustrate these properties in Fig. |aT] which shows the velocity er¬ 
ror from both schemes using either cubic spline or double hump cubic spline 
kernels. The performance of both schemes is similar, and as expected equa¬ 
tion |A1| is less sensitive to the choice of kernel, although it is still more 
accurate for the standard cubic-spline kernel. Since we have investigated 
the drag forces in the test-particle limit, in which exact angular momentum 
conservation is not an issue, the comparable accuracy justifies our choice of 
force law, which was essentially motivated by facilitating a direct compari¬ 
son to earlier work. However, we suggest equation [A^ should be used when 
the force from the dust on the gas cannot be neglected. 
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